#pragma once
#include <vector>
#include <algorithm>
using namespace std;
#include "point.h"
#include <Eigen/Sparse>
#include <Eigen/SparseLU>
#include <Eigen/PardisoSupport>
using namespace Eigen;

enum BASIS_TYPE//basis function type
{
	TWO_DIM_LINEAR,//2D linear
	TWO_DIM_QUADRATIC,//2D quadratic
	TWO_DIM_BILINEAR,//2D bilinear
	TWO_DIM_BIQUADRATIC//2D biquadratic
};

enum BOUNDARY_TYPE//boundary type
{
	Dirichlet,
	Neumann,
	Robin,
	Other
};



enum ELE_TYPE {//finite element type
	ELE_3 = 3,//Linear triangle element with 3 nodes
	ELE_4 = 4,//Bilinear element with 4 nodes
	ELE_6 = 6,//quadratic triangle element with 6 nodes
	ELE_9 = 9//Biquadratic element with 9 nodes
};

enum SOLVE_METHOD {
	DIREC_SparseLU,
	DIREC_SimplicialCholesky,
	DIREC_PardisoLU,
	DIREC_PardisoLDLT,
	ITERA_BiCGSTAB,

};

//One node for boundary
struct boundary_node
{
	BOUNDARY_TYPE type;//define the boundary type
	int glo_node_index;//define the global index of the boundary node
	vector<double(*)(CPoint, double t)> value;//define the boundary value, it is a vector because the unknown of one node may be more than one
	CPoint outer_normal_vector;//outer normal vector of the boundary used by Robin and Neumann conditions
	vector<int> mark;//define the which one should be used in the "value".
};

//One edge for boundary
struct boundary_edge
{
	BOUNDARY_TYPE type;//type of the boundary edge
	int ele_index;// the index of the element which contains the boundary edge
	pair<int,int> p;//the global node of the boundary edge
	vector<double(*)(CPoint,double t)> value;//define the boundary value of Neumann or Robin, it is defined by a function
	vector<double(*)(CPoint,double t)> r;//define the coefficient before the Dirichlet condition

	boundary_edge(int e, pair<int,int>p0) { ele_index = e, p = p0; }//construct function for boundary_edge
	boundary_edge() {};
	vector<int> mark;//define the which one should be used in the "value" and "r".
};

typedef Eigen::SparseMatrix<double> SpMat;//declares a sparse matrix type of double
typedef Eigen::VectorXd Rsv;//declares a right side value type
typedef Eigen::Triplet<double> Tri;//declares a triplet to stores the matrix index and values, then convert it to sparse matrix(with high performance)
typedef Eigen::PardisoLU   <SpMat > PLUSolver;
typedef Eigen::PardisoLDLT <SpMat> PLDLTSolver;

template<typename Scalar, typename StorageIndex>
SparseMatrix<Scalar, ColMajor, StorageIndex> sparse_stack_v(
	const SparseMatrix<Scalar, ColMajor, StorageIndex>& top,
	const SparseMatrix<Scalar, ColMajor, StorageIndex>& bottom)
{
	assert(top.cols() == bottom.cols());
	SparseMatrix<Scalar, ColMajor, StorageIndex> stacked;
	stacked.resize(top.rows() + bottom.rows(), top.cols());
	stacked.resizeNonZeros(top.nonZeros() + bottom.nonZeros());

	StorageIndex i = 0;

	for (StorageIndex col = 0; col < top.cols(); col++)
	{
		stacked.outerIndexPtr()[col] = i;

		for (StorageIndex j = top.outerIndexPtr()[col]; j < top.outerIndexPtr()[col + 1]; j++, i++)
		{
			stacked.innerIndexPtr()[i] = top.innerIndexPtr()[j];
			stacked.valuePtr()[i] = top.valuePtr()[j];
		}

		for (StorageIndex j = bottom.outerIndexPtr()[col]; j < bottom.outerIndexPtr()[col + 1]; j++, i++)
		{
			stacked.innerIndexPtr()[i] = (StorageIndex)top.rows() + bottom.innerIndexPtr()[j];
			stacked.valuePtr()[i] = bottom.valuePtr()[j];
		}
	}
	stacked.outerIndexPtr()[top.cols()] = i;
}

template<typename Scalar, typename StorageIndex>
SparseMatrix<Scalar, ColMajor, StorageIndex> sparse_stack_h(
	const SparseMatrix<Scalar, ColMajor, StorageIndex>& left,
	const SparseMatrix<Scalar, ColMajor, StorageIndex>& right)
{
	assert(left.rows() == right.rows());

	SparseMatrix<Scalar, ColMajor, StorageIndex> stacked;

	stacked.resize(left.rows(), left.cols() + right.cols());
	stacked.resizeNonZeros(left.nonZeros() + right.nonZeros());

	std::copy(left.innerIndexPtr(), left.innerIndexPtr() + left.nonZeros(), stacked.innerIndexPtr());
	std::copy(right.innerIndexPtr(), right.innerIndexPtr() + right.nonZeros(), stacked.innerIndexPtr() + left.nonZeros());

	std::copy(left.valuePtr(), left.valuePtr() + left.nonZeros(), stacked.valuePtr());
	std::copy(right.valuePtr(), right.valuePtr() + right.nonZeros(), stacked.valuePtr() + left.nonZeros());

	std::copy(left.outerIndexPtr(), left.outerIndexPtr() + left.cols(), stacked.outerIndexPtr());//dont need the last entry of A.outerIndexPtr() -- total length is AB.cols() + 1 = A.cols() + B.cols() + 1
	std::transform(right.outerIndexPtr(), right.outerIndexPtr() + right.cols() + 1, stacked.outerIndexPtr() + left.cols(), [&](StorageIndex i) { return i + left.nonZeros(); });
}

template<typename Scalar, typename StorageIndex>
void sparse_stack_h_inplace(
	SparseMatrix<Scalar, ColMajor, StorageIndex>& left,
	const SparseMatrix<Scalar, ColMajor, StorageIndex>& right)
{
	assert(left.rows() == right.rows());

	const StorageIndex leftcol = (StorageIndex)left.cols();
	const StorageIndex leftnz = (StorageIndex)left.nonZeros();

	left.conservativeResize(left.rows(), left.cols() + right.cols());
	left.resizeNonZeros(left.nonZeros() + right.nonZeros());

	std::copy(right.innerIndexPtr(), right.innerIndexPtr() + right.nonZeros(), left.innerIndexPtr() + leftnz);
	std::copy(right.valuePtr(), right.valuePtr() + right.nonZeros(), left.valuePtr() + leftnz);
	std::transform(right.outerIndexPtr(), right.outerIndexPtr() + right.cols() + 1, left.outerIndexPtr() + leftcol, [&](StorageIndex i) { return i + leftnz; });
}